Chapter 6 Diversity analysis

6.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

6.1.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.3 Antibiotics samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="2_Antibiotics") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.4 Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="3_Transplant1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.5 Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="4_Transplant2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.6 Post-Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.7 Post-Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

6.3 Permanovas

6.3.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

6.3.1 1. Are the wild populations similar?

6.3.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")

6.3.1.2 Number of samples used

[1] 27

beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts, q=1, dist=dist)

6.3.1.3 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000012 0.000012 0.0012    999  0.979
Residuals 25 0.257281 0.010291                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.977
Podarcis_muralis            0.97302                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.542719	0.2095041	6.625717	0.001
Residual	25	5.820951	0.7904959	NA	NA
Total	26	7.363669	1.0000000	NA	NA

6.3.1.4 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000048 0.0000476 0.0044    999  0.937
Residuals 25 0.270114 0.0108046                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.935
Podarcis_muralis            0.94763                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.918266	0.2608511	8.822682	0.001
Residual	25	5.435610	0.7391489	NA	NA
Total	26	7.353876	1.0000000	NA	NA

6.3.1.5 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03585 0.035847 2.4912    999  0.125
Residuals 25 0.35973 0.014389                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.125
Podarcis_muralis            0.12705                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.3218613	0.2162815	6.899207	0.001
Residual	25	1.1662981	0.7837185	NA	NA
Total	26	1.4881594	1.0000000	NA	NA

6.3.1.6 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.018367 0.018367 1.5597    999  0.246
Residuals 25 0.294402 0.011776                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                  Podarcis_liolepis Podarcis_muralis
Podarcis_liolepis                              0.234
Podarcis_muralis            0.22328                 

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0858578	0.172879	5.225323	0.051
Residual	25	0.4107775	0.827121	NA	NA
Total	26	0.4966352	1.000000	NA	NA

beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

6.3.2 2. Effect of acclimation

accli <- meta %>%
  filter(time_point == "1_Acclimation")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(accli))]
identical(sort(colnames(accli.counts)), sort(as.character(rownames(accli))))

accli_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation")

6.3.2.1 Number of samples used

[1] 27

beta_div_richness_accli<-hillpair(data=accli.counts, q=0)
beta_div_neutral_accli<-hillpair(data=accli.counts, q=1)
beta_div_phylo_accli<-hillpair(data=accli.counts, q=1, tree=genome_tree)
beta_div_func_accli<-hillpair(data=accli.counts, q=1, dist=dist)

6.3.2.2 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.11796 0.117959 12.963    999  0.003 **
Residuals 25 0.22748 0.009099                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.003
Hot_dry  0.0013711        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	1.639807	0.179834	5.481634	0.001
Residual	25	7.478640	0.820166	NA	NA
Total	26	9.118447	1.000000	NA	NA

6.3.2.3 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.07844 0.078443 5.2384    999  0.027 *
Residuals 25 0.37437 0.014975                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.025
Hot_dry  0.030815        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	1.947003	0.2306127	7.493387	0.001
Residual	25	6.495736	0.7693873	NA	NA
Total	26	8.442739	1.0000000	NA	NA

6.3.2.4 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.06739 0.067395 2.9532    999  0.086 .
Residuals 25 0.57052 0.022821                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.093
Hot_dry  0.098068        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	0.2441653	0.1224638	3.488854	0.019
Residual	25	1.7496100	0.8775362	NA	NA
Total	26	1.9937754	1.0000000	NA	NA

6.3.2.5 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02496 0.024955 0.6729    999   0.41
Residuals 25 0.92714 0.037085                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.39
Hot_dry   0.41979        

	Df	SumOfSqs	R2	F	Pr(>F)
Population	1	0.0279454	0.0248037	0.6358634	0.447
Residual	25	1.0987171	0.9751963	NA	NA
Total	26	1.1266624	1.0000000	NA	NA

beta_q0n_nmds_accli <- beta_div_richness_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_accli <- beta_div_neutral_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_accli <- beta_div_phylo_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_accli <- beta_div_func_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

6.3.3 3. Comparison between Wild and Acclimation

accli1 <- meta  %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(accli1))]
identical(sort(colnames(accli1.counts)),sort(as.character(rownames(accli1))))

accli1_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

6.3.3.1 Number of samples used

[1] 54

beta_div_richness_accli1<-hillpair(data=accli1.counts, q=0)
beta_div_neutral_accli1<-hillpair(data=accli1.counts, q=1)
beta_div_phylo_accli1<-hillpair(data=accli1.counts, q=1, tree=genome_tree)
beta_div_func_accli1<-hillpair(data=accli1.counts, q=1, dist=dist)

6.3.3.1.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.05014 0.050145 6.2252    999  0.014 *
Residuals 52 0.41886 0.008055                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                0_Wild 1_Acclimation
0_Wild                         0.012
1_Acclimation 0.015808              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.6172653	0.0360987	3.933397	0.001
species	1	2.8279677	0.1653842	18.020647	0.001
individual	25	9.5739861	0.5599025	2.440331	0.001
Residual	26	4.0801621	0.2386146	NA	NA
Total	53	17.0993812	1.0000000	NA	NA

6.3.3.1.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0199 0.0199035 2.1213    999  0.157
Residuals 52 0.4879 0.0093827                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.155
1_Acclimation 0.15128              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.9050519	0.0541893	6.651487	0.001
species	1	3.3236300	0.1989999	24.426315	0.001
individual	25	8.9352276	0.5349902	2.626702	0.001
Residual	26	3.5377576	0.2118206	NA	NA
Total	53	16.7016671	1.0000000	NA	NA

6.3.3.1.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01334 0.013340 0.6524    999  0.448
Residuals 52 1.06332 0.020449                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.453
1_Acclimation 0.42294              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.2890434	0.0766494	7.532050	0.001
species	1	0.3508889	0.0930498	9.143655	0.001
individual	25	2.1332925	0.5657133	2.223620	0.004
Residual	26	0.9977533	0.2645874	NA	NA
Total	53	3.7709782	1.0000000	NA	NA

6.3.3.1.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0123 0.012300 0.4817    999   0.53
Residuals 52 1.3277 0.025533                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.525
1_Acclimation 0.49073              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	0.0448774	0.0269021	2.355512	0.172
species	1	0.0973005	0.0583275	5.107077	0.046
individual	25	1.0306426	0.6178264	2.163841	0.071
Residual	26	0.4953546	0.2969440	NA	NA
Total	53	1.6681751	1.0000000	NA	NA

beta_richness_nmds_accli1 <- beta_div_richness_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_accli1 <- beta_div_neutral_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_accli1 <- beta_div_phylo_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_accli1 <- beta_div_func_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

6.3.4 4. Do the antibiotics work?

6.3.4.1 Acclimation vs antibiotics

treat <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))

treat_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

6.3.4.2 Number of samples used

[1] 50

beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)

6.3.4.2.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq     F N.Perm Pr(>F)  
Groups     1 0.025318 0.0253178 6.021    999  0.014 *
Residuals 48 0.201837 0.0042049                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.018
2_Antibiotics      0.017817              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.888584	0.0949462	6.361098	0.001
species	1	2.117109	0.1064350	7.130814	0.001
individual	25	9.353701	0.4702455	1.260199	0.003
Residual	22	6.531709	0.3283734	NA	NA
Total	49	19.891103	1.0000000	NA	NA

6.3.4.2.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.039587 0.039587 6.8387    999  0.009 **
Residuals 48 0.277854 0.005789                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                        0.01
2_Antibiotics      0.011886              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	2.024181	0.1063620	9.051981	0.001
species	1	2.853103	0.1499183	12.758858	0.001
individual	25	9.234189	0.4852168	1.651783	0.001
Residual	22	4.919584	0.2585029	NA	NA
Total	49	19.031057	1.0000000	NA	NA

6.3.4.2.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.58372 0.58372 35.413    999  0.001 ***
Residuals 48 0.79119 0.01648                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                       0.001
2_Antibiotics    2.9795e-07              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.8065206	0.2113909	18.636551	0.001
species	1	0.7903334	0.0924813	8.153292	0.001
individual	25	3.8164689	0.4465860	1.574869	0.007
Residual	22	2.1325541	0.2495419	NA	NA
Total	49	8.5458771	1.0000000	NA	NA

6.3.4.2.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.18591 0.185914 5.0679    999  0.027 *
Residuals 48 1.76088 0.036685                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 2_Antibiotics
1_Acclimation                        0.03
2_Antibiotics      0.028989              

	Df	SumOfSqs	R2	F	Pr(>F)
time_point	1	1.8020952	0.3750193	33.6195614	0.001
species	1	0.0031247	0.0006503	0.0582938	0.831
individual	25	1.8208629	0.3789249	1.3587875	0.192
Residual	22	1.1792568	0.2454055	NA	NA
Total	49	4.8053396	1.0000000	NA	NA

beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat <- beta_div_func_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

6.3.5 5. Does the FMT work?

6.3.5.1 Comparison between FMT2 vs Post-FMT2

#Create newID to identify duplicated samples
transplants_metadata<-sample_metadata%>%
  mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)

transplant3<-transplants_metadata%>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")%>%
  column_to_rownames("newID")

transplant3_nmds <- transplants_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

full_counts<-temp_genome_counts %>%
    t()%>%
    as.data.frame()%>%
    rownames_to_column("Tube_code")%>%
    full_join(transplants_metadata,by = join_by(Tube_code == Tube_code))

transplant3_counts<-full_counts %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2") %>%
  subset(select=-c(315:324)) %>%
  column_to_rownames("newID")%>%
  subset(select=-c(1))%>%
  t() %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)

identical(sort(colnames(transplant3_counts)),sort(as.character(rownames(transplant3))))

6.3.5.2 Number of samples used

[1] 49

beta_div_richness_transplant3<-hillpair(data=transplant3_counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3_counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3_counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3_counts, q=1, dist=dist)

#Arrange of metadata dataframe
transplant3_arrange<-transplant3[labels(beta_div_neutral_transplant3$S),]

6.3.5.2.1 Richness

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.180473	0.0855095	6.984555	0.001
time_point	1	0.860906	0.0623612	5.093759	0.001
type	1	1.459433	0.1057165	8.635089	0.001
individual	24	6.755100	0.4893170	1.665341	0.001
Residual	21	3.549250	0.2570959	NA	NA
Total	48	13.805162	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 1.4169018 5.739828 0.15622903   0.001      0.003   *
2   Control vs Hot_control  1 2.0940966 8.509112 0.21005427   0.001      0.003   *
3 Treatment vs Hot_control  1 0.3004618 1.265034 0.04179854   0.144      0.432    

6.3.5.2.2 Neutral

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	1.2800927	0.0939787	8.796453	0.001
time_point	1	0.9350566	0.0686477	6.425458	0.001
type	1	1.9135997	0.1404879	13.149743	0.001
individual	24	6.4363516	0.4725281	1.842870	0.001
Residual	21	3.0559984	0.2243577	NA	NA
Total	48	13.6210990	1.0000000	NA	NA

                     pairs Df SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 1.8758788  8.282671 0.21084796   0.001      0.003   *
2   Control vs Hot_control  1 2.4396317 10.635546 0.24945256   0.001      0.003   *
3 Treatment vs Hot_control  1 0.3158428  1.394345 0.04587515   0.115      0.345    

6.3.5.2.3 Phylogenetic

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.1400466	0.0952654	6.956436	0.001
time_point	1	0.1138047	0.0774145	5.652939	0.001
type	1	0.1432667	0.0974558	7.116383	0.001
individual	24	0.6501795	0.4422784	1.345663	0.046
Residual	21	0.4227709	0.2875859	NA	NA
Total	48	1.4700683	1.0000000	NA	NA

                     pairs Df  SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.14387705 5.735321 0.15612552   0.001      0.003   *
2   Control vs Hot_control  1 0.22715701 9.044894 0.22036587   0.001      0.003   *
3 Treatment vs Hot_control  1 0.04648319 1.704277 0.05550617   0.123      0.369    

6.3.5.2.4 Functional

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0092808	0.0077189	0.4182529	0.528
time_point	1	-0.0061674	-0.0051295	-0.2779456	0.883
type	1	0.0831052	0.0691191	3.7452726	0.101
individual	24	0.6501528	0.5407359	1.2208414	0.349
Residual	21	0.4659767	0.3875556	NA	NA
Total	48	1.2023481	1.0000000	NA	NA

                     pairs Df    SumsOfSqs     F.Model           R2 p.value p.adjusted sig
1     Control vs Treatment  1  0.078539743  4.59293783  0.129040706   0.089      0.267    
2   Control vs Hot_control  1  0.052468954  2.13675422  0.062593948   0.180      0.540    
3 Treatment vs Hot_control  1 -0.002340352 -0.07432315 -0.002569452   0.857      1.000    

beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == newID))

p0<-beta_richness_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylo_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_func_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

6.3.5.3 Comparison between the different experimental time points

6.3.5.3.1 Comparison againts Wild samples

transplants_metadata<-sample_metadata%>%
  mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)

transplant4<-transplants_metadata%>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2" |time_point == "0_Wild")%>%
  column_to_rownames("newID")

transplant4_pairwise <- transplants_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2"|time_point == "0_Wild")


full_counts<-temp_genome_counts %>%
  t()%>%
  as.data.frame()%>%
  rownames_to_column("Tube_code")%>%
  full_join(transplants_metadata,by = join_by(Tube_code == Tube_code))

transplant4_counts<-full_counts %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2"|time_point == "0_Wild") %>%
  subset(select=-c(315:324)) %>%
  column_to_rownames("newID")%>%
  subset(select=-c(1))%>%
  t() %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)

identical(sort(colnames(transplant4_counts)),sort(as.character(rownames(transplant4))))

beta_div_neutral_transplant4<-hillpair(data=transplant4_counts, q=1)

The estimated time for calculating the 2850 pairwise combinations is 13 seconds.

transplant4_beta_div_neutral<-as.matrix(beta_div_neutral_transplant4$S) #convert the pairwise comparisons into a matrix

transplant4_beta_div_neutral_df <- as.data.frame(as.table(as.matrix(transplant4_beta_div_neutral))) #convert the matrix to a data frame with pairwise values

transplant4_beta_div_neutral_df <- transplant4_beta_div_neutral_df %>% #remove the duplicated pairwise comparisons
  filter(Var1 != Var2) %>%
  mutate(pair = pmap_chr(list(Var1, Var2), ~paste(sort(c(..1, ..2)), collapse = "_"))) %>%
  distinct(pair, .keep_all = TRUE) %>%
  select(Sample_1 = Var1, Sample_2 = Var2, Distance = Freq) %>%
  arrange(Sample_1, Sample_2)


transplant4_beta_div_neutral_df <- transplant4_beta_div_neutral_df %>% #keep only the pairwise comparisons between the same individual
  mutate(Sample_1_part = sub("^[^_]*_", "", Sample_1),
         Sample_2_part = sub("^[^_]*_", "", Sample_2)) %>%
  filter(Sample_1_part == Sample_2_part) 
#%>%
  #select(Sample_1, Sample_2, Distance)


transplant4_beta_div_neutral_met<-transplant4_beta_div_neutral_df %>% #merge with the metadata associated to the pairwise comparisons for Sample_1
  inner_join(transplant4_pairwise, by = join_by(Sample_1 == "newID"))

# Extract the part before the first `_` for Sample_1 and Sample_2
transplant4_beta_div_neutral_met$Sample_1_part <- sub("_.*", "", transplant4_beta_div_neutral_met$Sample_1)

# Create a named vector for matching Tube_code with Time_point
time_point_map <- setNames(transplant4_beta_div_neutral_met$time_point, transplant4_beta_div_neutral_met$Tube_code)

# Function to replace the part before the first _ with the corresponding Time_point
replace_with_time_point <- function(sample_part, tube_code, time_point_map) {
  if (tube_code %in% names(time_point_map)) {
    return(time_point_map[tube_code])
  } else {
    return(sample_part)
  }
}

# Apply the function to Sample_1_part
transplant4_beta_div_neutral_met$Sample_1_part <- mapply(replace_with_time_point,
                                                         transplant4_beta_div_neutral_met$Sample_1_part,
                                                         transplant4_beta_div_neutral_met$Tube_code,
                                                         list(time_point_map))

transplant4_beta_div_neutral_met<-transplant4_beta_div_neutral_met %>%
  subset(select=-c(6:16))

transplant4_beta_div_neutral_met$Sample_2_part <- sub("_.*", "", transplant4_beta_div_neutral_met$Sample_2)


transplant4_beta_div_neutral_met<-transplant4_beta_div_neutral_met %>% #merge with the metadata associated to the pairwise comparisons for Sample_2
  inner_join(transplant4_pairwise, by = join_by(Sample_2 == "newID"))

# Create a named vector again for matching Tube_code with Time_point
time_point_map <- setNames(transplant4_beta_div_neutral_met$time_point, transplant4_beta_div_neutral_met$Tube_code)


replace_with_time_point <- function(sample_part, tube_code, time_point_map) {
  if (tube_code %in% names(time_point_map)) {
    return(time_point_map[tube_code])
  } else {
    return(sample_part)
  }
}

# Apply the function to Sample_2_part
transplant4_beta_div_neutral_met$Sample_2_part <- mapply(replace_with_time_point,
                                                    transplant4_beta_div_neutral_met$Sample_2_part,
                                                         transplant4_beta_div_neutral_met$Tube_code,
                                                         list(time_point_map))


# Create a new variable to plot the distances between the time_points
transplant4_beta_div_neutral_met$comparison <- paste(transplant4_beta_div_neutral_met$Sample_1_part, "-",transplant4_beta_div_neutral_met$Sample_2_part)

##plot the differences between wild and post-transplant and transplant and post-transplant
transplant4_beta_div_neutral_met  %>% 
  ggplot(aes(x=type, y=Distance, color=type, fill=type))+
    geom_boxplot() +
    scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
    scale_fill_manual(name="Type",
                      breaks=c("Control", "Hot_control", "Treatment"),
                      labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                      values=c("#4477AA50","#d57d2c50","#76b18350")) +
    scale_x_discrete(labels=c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
    theme_minimal() +
    facet_wrap(~comparison, labeller = labeller(comparison = c("0_Wild - 6_Post-FMT2" = "Wild-Treatment", "0_Wild - 4_Transplant2" = "Wild-Transplant", "6_Post-FMT2 - 4_Transplant2"="Treatment-Transplant"))) +
    theme(
    axis.text.x = element_text(angle = 45, hjust = 1))

6.3.5.3.2 Comparison againts Acclimation samples

#Create newID to identify duplicated samples
transplants_metadata<-sample_metadata%>%
  mutate(Tube_code=str_remove_all(Tube_code, "_a"))
transplants_metadata$newID <- paste(transplants_metadata$Tube_code, "_", transplants_metadata$individual)

transplant5<-transplants_metadata%>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2" |time_point == "1_Acclimation")%>%
  column_to_rownames("newID")

transplant5_pairwise <- transplants_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2"|time_point == "1_Acclimation")


full_counts_new<-temp_genome_counts %>%
  t()%>%
  as.data.frame()%>%
  rownames_to_column("Tube_code")%>%
  full_join(transplants_metadata,by = join_by(Tube_code == Tube_code))

transplant5_counts<-full_counts %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2"|time_point == "1_Acclimation") %>%
  subset(select=-c(315:324)) %>%
  column_to_rownames("newID")%>%
  subset(select=-c(1))%>%
  t() %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)

identical(sort(colnames(transplant5_counts)),sort(as.character(rownames(transplant5))))


beta_div_neutral_transplant5<-hillpair(data=transplant5_counts, q=1)

The estimated time for calculating the 2850 pairwise combinations is 13 seconds.

transplant5_beta_div_neutral<-as.matrix(beta_div_neutral_transplant5$S) #convert the pairwise comparisons into a matrix

transplant5_beta_div_neutral_df <- as.data.frame(as.table(as.matrix(transplant5_beta_div_neutral))) #convert the matrix to a data frame with pairwise values

transplant5_beta_div_neutral_df <- transplant5_beta_div_neutral_df %>% #remove the duplicated pairwise comparisons
  filter(Var1 != Var2) %>%
  mutate(pair = pmap_chr(list(Var1, Var2), ~paste(sort(c(..1, ..2)), collapse = "_"))) %>%
  distinct(pair, .keep_all = TRUE) %>%
  select(Sample_1 = Var1, Sample_2 = Var2, Distance = Freq) %>%
  arrange(Sample_1, Sample_2)


transplant5_beta_div_neutral_df <- transplant5_beta_div_neutral_df %>% #keep only the pairwise comparisons between the same individual
  mutate(Sample_1_part = sub("^[^_]*_", "", Sample_1),
         Sample_2_part = sub("^[^_]*_", "", Sample_2)) %>%
  filter(Sample_1_part == Sample_2_part) 
#%>%
  #select(Sample_1, Sample_2, Distance)


transplant5_beta_div_neutral_met<-transplant5_beta_div_neutral_df %>% #merge with the metadata associated to the pairwise comparisons for Sample_1
  inner_join(transplant5_pairwise, by = join_by(Sample_1 == "newID"))

# Extract the part before the first `_` for Sample_1 and Sample_2
transplant5_beta_div_neutral_met$Sample_1_part <- sub("_.*", "", transplant5_beta_div_neutral_met$Sample_1)

# Create a named vector for matching Tube_code with Time_point
time_point_map_new <- setNames(transplant5_beta_div_neutral_met$time_point, transplant5_beta_div_neutral_met$Tube_code)

# Function to replace the part before the first _ with the corresponding Time_point
replace_with_time_point_new <- function(sample_part, tube_code, time_point_map_new) {
  if (tube_code %in% names(time_point_map_new)) {
    return(time_point_map_new[tube_code])
  } else {
    return(sample_part)
  }
}

# Apply the function to Sample_1_part
transplant5_beta_div_neutral_met$Sample_1_part <- mapply(replace_with_time_point_new,
                                                         transplant5_beta_div_neutral_met$Sample_1_part,
                                                         transplant5_beta_div_neutral_met$Tube_code,
                                                         list(time_point_map_new))

transplant5_beta_div_neutral_met<-transplant5_beta_div_neutral_met %>%
  subset(select=-c(6:16))

transplant5_beta_div_neutral_met$Sample_2_part <- sub("_.*", "", transplant5_beta_div_neutral_met$Sample_2)


transplant5_beta_div_neutral_met<-transplant5_beta_div_neutral_met %>% #merge with the metadata associated to the pairwise comparisons for Sample_2
  inner_join(transplant5_pairwise, by = join_by(Sample_2 == "newID"))

# Create a named vector again for matching Tube_code with Time_point
time_point_map_new <- setNames(transplant5_beta_div_neutral_met$time_point, transplant5_beta_div_neutral_met$Tube_code)


replace_with_time_point_new <- function(sample_part, tube_code, time_point_map_new) {
  if (tube_code %in% names(time_point_map_new)) {
    return(time_point_map_new[tube_code])
  } else {
    return(sample_part)
  }
}

# Apply the function to Sample_2_part
transplant5_beta_div_neutral_met$Sample_2_part <- mapply(replace_with_time_point_new,
                                                         transplant5_beta_div_neutral_met$Sample_2_part,
                                                         transplant5_beta_div_neutral_met$Tube_code,
                                                         list(time_point_map_new))


# Create a new variable to plot the distances between the time_points
transplant5_beta_div_neutral_met$comparison <- paste(transplant5_beta_div_neutral_met$Sample_1_part, "-",transplant5_beta_div_neutral_met$Sample_2_part)


##plot the differences between acclimation and post-transplant and transplant and post-transplant

# Reorder the 'comparison' factor
transplant5_beta_div_neutral_met$comparison <- factor(
  transplant5_beta_div_neutral_met$comparison,
  levels = c("4_Transplant2 - 1_Acclimation",  # Acclimation-Transplant first
             "6_Post-FMT2 - 4_Transplant2",    # Transplant-Treatment second
             "6_Post-FMT2 - 1_Acclimation")    # Treatment-Acclimation last
)

# Create the plot with the specified facet order
transplant5_beta_div_neutral_met %>% 
  ggplot(aes(x=type, y=Distance, color=type, fill=type)) +
  geom_boxplot() +
  scale_color_manual(name="Type",
                     breaks=c("Control", "Hot_control", "Treatment"),
                     labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                     values=c("#4477AA","#d57d2c","#76b183")) +
  scale_fill_manual(name="Type",
                    breaks=c("Control", "Hot_control", "Treatment"),
                    labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                    values=c("#4477AA50","#d57d2c50","#76b18350")) +
  scale_x_discrete(labels=c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
  theme_minimal() +
  facet_wrap(~comparison, 
             labeller = labeller(comparison = c(
               "4_Transplant2 - 1_Acclimation" = "Acclimation-Transplant", 
               "6_Post-FMT2 - 4_Transplant2" = "Transplant-Treatment", 
               "6_Post-FMT2 - 1_Acclimation" = "Treatment-Acclimation"
             )),
             nrow = 1, ncol = 3) +  # Single row with 3 columns
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  ylab("Beta diversity dissimilarity distance") +
  xlab("Type")

6.3.6 7. Are there differences between the control and the treatment group?

6.3.6.1 After 1 week –> Post-FMT1

post1 <- meta %>%
  filter(time_point == "5_Post-FMT1")

post1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post1))]
identical(sort(colnames(post1.counts)),sort(as.character(rownames(post1))))

post1_nmds <- sample_metadata %>%
  filter(time_point == "5_Post-FMT1")

6.3.6.2 Number of samples used

[1] 26

beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post1_arrange<-post1[labels(beta_div_neutral_post1$S),]

6.3.6.2.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.017675 0.0088373 2.3825    999  0.095 .
Residuals 23 0.085312 0.0037092                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0030000     0.647
Hot_control 0.0068795                 0.215
Treatment   0.6248469   0.2084296          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.6340254	0.0768024	2.065607	0.004
type	1	0.5615418	0.0680222	1.829461	0.009
Residual	23	7.0597099	0.8551754	NA	NA
Total	25	8.2552771	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.5615418 1.729004 0.1033537   0.020      0.060    
2   Control vs Hot_control  1 0.8438429 2.793772 0.1486541   0.001      0.003   *
3 Treatment vs Hot_control  1 0.3734921 1.268929 0.0779971   0.103      0.309    

6.3.6.2.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.011001 0.0055005 0.6303    999   0.54
Residuals 23 0.200714 0.0087267                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.21300     0.961
Hot_control 0.21166                 0.448
Treatment   0.95468     0.43604          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.7907904	0.1076445	3.056657	0.001
type	1	0.6051778	0.0823784	2.339205	0.010
Residual	23	5.9503501	0.8099772	NA	NA
Total	25	7.3463184	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.6051778 2.250849 0.13047758   0.017      0.051    
2   Control vs Hot_control  1 1.0528902 4.143637 0.20570451   0.001      0.003   *
3 Treatment vs Hot_control  1 0.4150076 1.637268 0.09840968   0.044      0.132    

6.3.6.2.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00440 0.0021994 0.1369    999  0.926
Residuals 23 0.36941 0.0160614                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.94500     0.719
Hot_control 0.91505                 0.775
Treatment   0.63312     0.73046          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0560850	0.0531376	1.3149967	0.262
type	1	0.0184254	0.0174571	0.4320099	0.791
Residual	23	0.9809570	0.9294053	NA	NA
Total	25	1.0554673	1.0000000	NA	NA

                     pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.01842535 0.4144162 0.02688498   0.782      1.000    
2   Control vs Hot_control  1 0.05987967 1.7387847 0.09802164   0.106      0.318    
3 Treatment vs Hot_control  1 0.03212966 0.6477782 0.04139746   0.690      1.000    

6.3.6.2.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.00400 0.0020014 0.145    999  0.873
Residuals 23 0.31753 0.0138057                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.61200     0.745
Hot_control 0.59817                 0.849
Treatment   0.75141     0.83718          

	Df	SumOfSqs	R2	F	Pr(>F)
species	1	0.0024979	0.0033024	0.0900845	0.625
type	1	0.1161466	0.1535542	4.1887855	0.076
Residual	23	0.6377435	0.8431434	NA	NA
Total	25	0.7563879	1.0000000	NA	NA

                     pairs Df  SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.11614656 4.724791 0.23953568   0.073      0.219    
2   Control vs Hot_control  1 0.05000930 1.704826 0.09629160   0.229      0.687    
3 Treatment vs Hot_control  1 0.01235859 0.423812 0.02747777   0.471      1.000    

beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.6.3 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

6.3.6.4 Number of samples used

[1] 27

beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post2_arrange<-post2[labels(beta_div_neutral_post2$S),]

6.3.6.4.1 Richness

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.002011 0.0010056 0.1982    999  0.818
Residuals 24 0.121775 0.0050740                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.71700     0.807
Hot_control 0.67789                 0.571
Treatment   0.79246     0.59820          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	1.504341	0.1967776	2.939822	0.001
Residual	24	6.140538	0.8032224	NA	NA
Total	26	7.644879	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Treatment vs Control  1 0.6463814 2.560441 0.1379515   0.001      0.003   *
2 Treatment vs Hot_control  1 0.4796256 1.916520 0.1069694   0.003      0.009   *
3   Control vs Hot_control  1 1.1305044 4.268317 0.2105906   0.001      0.003   *

6.3.6.4.2 Neutral

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008262 0.0041311 0.8024    999  0.479
Residuals 24 0.123559 0.0051483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.42400     0.682
Hot_control 0.44675                 0.256
Treatment   0.65989     0.25095          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	1.923807	0.2603795	4.224537	0.001
Residual	24	5.464666	0.7396205	NA	NA
Total	26	7.388473	1.0000000	NA	NA

                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Treatment vs Control  1 1.0227481 4.648335 0.2251191   0.001      0.003   *
2 Treatment vs Hot_control  1 0.5010202 2.206532 0.1211945   0.002      0.006   *
3   Control vs Hot_control  1 1.3619424 5.771031 0.2650785   0.001      0.003   *

6.3.6.4.3 Phylogenetic

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000407 0.0002034 0.0487    999  0.969
Residuals 24 0.100305 0.0041794                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.95600     0.857
Hot_control 0.93765                 0.775
Treatment   0.83933     0.76015          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	0.1594363	0.2042241	3.079623	0.001
Residual	24	0.6212564	0.7957759	NA	NA
Total	26	0.7806927	1.0000000	NA	NA

                     pairs Df  SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Treatment vs Control  1 0.05927454 2.382025 0.1295845   0.033      0.099    
2 Treatment vs Hot_control  1 0.06906280 2.722460 0.1454115   0.004      0.012   .
3   Control vs Hot_control  1 0.11081709 4.043656 0.2017424   0.002      0.006   *

6.3.6.4.4 Functional

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01126 0.0056302 0.2861    999  0.811
Residuals 24 0.47233 0.0196806                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.53100     0.680
Hot_control 0.48255                 0.813
Treatment   0.60116     0.75643          

	Df	SumOfSqs	R2	F	Pr(>F)
type	2	-0.0038724	-0.0056213	-0.0670788	0.92
Residual	24	0.6927468	1.0056213	NA	NA
Total	26	0.6888744	1.0000000	NA	NA

                     pairs Df    SumsOfSqs     F.Model           R2 p.value p.adjusted sig
1     Treatment vs Control  1 -0.008527330 -0.46290555 -0.029793572   0.856          1    
2 Treatment vs Hot_control  1 -0.001648721 -0.04717131 -0.002956924   0.908          1    
3   Control vs Hot_control  1  0.004367477  0.13147026  0.008149924   0.679          1    

beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")